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Abstract 

The  objective  of  the  current  project  is  the  development  of  the  fundamentals  of  a  novel  two-scale 
multiscale  computational  method  for  the  nonlinear  damage  and  failure  analysis  of  3D  woven 
fiber  composites  under  ballistic  loading.  Since  material  behavior  is  detennined  by  its 
microstructure,  it  is  essential  to  accurately  model  the  physics  at  that  scale.  The  macroscale 
analysis  provides  a  useful  insight  into  the  underlying  high  strain  rate  physics  which  is  essential  in 
modeling  the  lower  micro-scale.  In  particular  a  rate  dependent  constitutive  approach  is  being 
developed  coupled  with  continuum  damage  mechanics  suitable  for  polymer  materials.  The  effect 
of  contact  parameters  on  the  underlying  damage  processes  is  being  studied  and  worked  on.  We 
further  develop  a  material  model  suitable  particularly  for  loading  of  composites  in  the  high 
strain  rate  regime.  This  is  a  significant  development  from  the  previous  model  where  strain 
rate  sensitivity  is  a-priori  postulated  for  the  matrix  dominated  modes  in  the  small  strain 
framework.  We  focused  on  developing  a  general  homogenized  anisotropic  material  model 
and  obtained  results  which  can  be  implemented  in  a  finite  element  framework  for  high  strain 
rate  loading. 
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1.  Introduction 

Ballistic  penetration  of  composites  involves  a  number  of  nonlinearities  acting  together  stemming 
from  high  strain  rates,  contact-impact  and  evolving  damage.  The  membrane  models  (Leigh 
Phoenix  and  Porwal  2003;  Nadler  et  al.  2006),  usually  developed  for  fabric  systems,  are 
inapplicable  in  the  case  of  3D  orthogonal  woven  composites  (3D-OWC)  because  of  elastic 
membrane  assumptions,  inability  to  incorporate  strain  rate  effects,  limitations  pertaining  to 
simple  stress/strain  based  failure  laws  and  lack  of  systematic  framework  to  complicated 
geometries  and  large  deformations.  This  limits  the  usability  of  these  models  for  advanced 
composites  like  3D-OWC.  On  the  other  hand,  a  microscale  simulation  with  resolution  of 
individual  fiber  filament  is  impractical  due  to  enormous  computational  resources  needed.  This 
has  led  to  various  homogenization  schemes  both  for  material  properties  and  damage  laws.  These 
methods  increasingly  popular  since  the  middle  of  the  last  decade  have  severe  technical 
difficulties.  Some  of  these  methods  try  to  drastically  cut  computing  time  by  unit  cell  level 
homogenization  (Baheieldin  et  al.  2004;  Bahei-El-Din  and  Zikry  2003;  Pankow  et  al.  2011; 
Pankow  et  al.  2011).  The  periodicity  of  the  unit  cell  makes  it  an  ideal  candidate  for  the  classical 
representative  volume  element  (RVE).  However,  this  assumption  is  invalid  for  short  wavelength 
regimes  such  as  a  pulse-loading  event  like  ballistic  impact.  Oskay  and  Fish  (2007)  sidestep  the 
problem  by  applying  asymptotic  homogenization  at  the  level  of  fiber  and  filament  through  an 
eigendefonnation  formulation.  Unfortunately  their  analysis  neglects  strain  rate  and  material 
nonlinearities  in  addition  to  requiring,  often  intractable,  mathematical  manipulations.  In  addition, 
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using  the  fiber  filament  and  surrounding  matrix  as  RVE  is  questionable  since  the  architecture  of 
an  actual  3D-OWC  is  far  more  complex.  It  must  also  be  noted  that  most  of  damage  laws  have  not 
been  compared  with  experimental  penetration  events. 

This  issue  is  recently  addressed  by  Sun  et  al.  (2009)  and  Jia  et  al.  (2010).  A  finite  element 
analysis  is  performed  with  a  critical  damage  area  (CD A)  theory  at  the  unit  cell  level  and  results 
are  compared  with  their  own  experiments.  Although  reasonable  match  between  experimental  and 
predicted  exit  velocities  were  reported,  many  of  the  material  and  damage  parameters  were  taken 
ad  hoc  without  a  physical  explanation.  We  would  like  to  point  out  that  developing  models  with  a 
large  number  of  parameters  fit  to  match  exit  velocities  only  is  not  a  viable  modeling  strategy  due 
to  its  sensitivity  to  mesh  distortion  parameters,  element  deletion  and  the  accuracy  of  the  contact- 
impact  conditions.  In  addition,  obtaining  exact  dynamic  material  parameters  for  every 
mechanical  phenomena  proposed  in  the  model  would  require  enormous  amount  of  testing  and 
calibration.  In  the  current  work  we  address  many  of  these  shortcomings.  We  incorporate  the 
nonlinearities  in  a  systematic  manner  employing  fewer  parameters  which  can  be  independently 
obtained  from  experimentation  of  constituents  of  the  composite.  In  addition,  we  employ  a 
damage  model  which  contains  fewer  parameters  whose  effects  on  the  constitutive  response  are 
easier  to  understand  and  calibrate,  are  more  intuitive  to  use  and  capture  the  physical 
consequences  of  the  various  damage  modes. 


Figure  1:  Spatial  hierarchichal  scales  in  a  typical  3D-OWC  applications. 

Constitutive  performance  of  heterogeneous  solids  cannot  be  accurately  predicted  if  the  effect 
of  microstructure  is  not  incorporated  into  the  constitutive  model.  Micro  structural  effects  are 
pronounced  in  composite  materials  because  of  the  presence  of  many  different  constituents  with 
widely  different  properties.  Such  a  varied  composition  is  useful  to  impart  desirable  properties 
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like  high  strength  to  weight  ratio  and  resistance  to  impact  damage.  However,  the  variation  of 
microstructure  exists  at  various  length  scales,  see  Figure  1. 

It  is  very  clear  that  computational  load  steadily  increases  as  we  move  up  spatial  resolution. 
Within  the  broad  division  of  atomic  and  continuum  scales  their  lie  a  number  of  length  scales.  We 
will  like  to  point  out  that  the  resolution  necessary  to  obtain  a  reasonably  accurate  response  is 
greatly  influenced  by  the  loading  of  the  structure.  As  the  impact  energy  increases,  the  resolution 
needed  to  accurately  capture  the  response  drops  sharply.  This  is  primarily  due  to  the  appearance 
of  distributed  micro  scale  damage  as  well  as  short  wavelength  waves  which  can  be  smaller  than 
the  RVE  of  the  medium.  Therefore,  it  is  inevitable  that  as  the  loading  becomes  more  severe,  for 
the  same  time  period  of  computation,  the  computing  requirement  greatly  increases  forcing  us  to 
make  an  appropriate  choice  of  the  length  scales.  We  characterize  the  scale  that  lies  between  the 
unit  cells  (cm)  to  the  fiber  filament  level  (am)  as  the  meso  scale.  We  use  this  scale  to  split  the 
sources  of  nonlinearities  -  inelasticity  and  damage  in  a  systematic  manner  which  is  both  simple 
as  well  as  deals  with  lesser  number  of  parameters.  We  exploit  two  simple  characteristics  of  3D- 
OWC  -  (1)  delamination  mode  of  damage  is  suppressed  for  localized  loading  due  to  3D  integral 
weave  in  one  layer  of  3D-OWC  and  (2)  closely  clustered  resin  impregnated  fibers  (henceforth 
we  call  it  yarns)  can  be  modeled  separately  from  the  resin  rich  regions  and  assumed  to  be  solely 
responsible  for  damage  accumulation  through  various  anisotropic  failure  modes  due  to  weak 
interfaces  and  stress  concentrators.  This  meso  level  scale  separation  reduces  computing  costs 
greatly,  obviates  the  theoretical  problems  with  conventional  multiscale  methods  and  helps  to 
systematically  deal  with  nonlinearities  according  to  the  physics  of  the  separated  regions. 

We  develop  a  partitioned  meso  scale  model  for  high  strain  rate  modeling  of  ballistic  impact 
on  three-dimensional  orthogonal  woven  composite.  In  this  model,  the  entire  material  nonlinearity 
stemming  from  large  deformation  and  strain  rate  effects  before  damage  initiation  was  borne  out 
by  the  polymeric  resin.  The  anisotropic  yams  are  modeled  using  extensions  of  quasi-static  small 
strain  models  with  strain  rate  sensitivity  empirically  appended.  We  further  extended  the  model  of 
the  yarn  into  a  large  defonnation  model. 

Modeling  of  large  strain  mechanics  has  been  attempted  before  for  anisotropic  composites 
(Simo  1987,  Matzenmiller  1995).  However,  they  have  not  been  extended  systematically  for  high 
strain  rate  loading  of  anisotropic  composites  where  important  material  parameter  like  Gruneisen 
tensor  and  damage  modes  becomes  dominant  (Anderson  et  al.  1990,  1994).  We  propose  a 
damage  law  based  on  consistent  thermomechanics  for  a  typical  unidirectional  composite,  which 
can  be  easily  extended  for  composites  with  greater  degree  of  anisotropy. 

The  report  is  organized  as  follows.  In  section  2  we  provide  the  details  of  the 
micromechanical  constitutive  model  of  resin.  In  section  3  constitutive  and  damage  model  for 
yam  is  outlined.  In  section  4  the  yam  model  is  further  extended  to  accommodate  finite 
deformation  along  with  a  yarn  level  damage  model  for  high  strain  rate  loading.  In  section  5  we 
discuss  the  simulation  results  followed  by  conclusions  in  section  6. 
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2.  Constitutive  modeling  of  resin 

The  polymer  resin  considered  in  the  current  model  is  epoxy  which  is  well  known  to  exhibit 
strong  rate  dependence  (Jordan  et  al.  2008),  not  considered  in  detail  in  any  previous  model  for 
ballistic  penetration  of  3D-OWC.  We  have  adopted  the  large  deformation,  high  strain  rate  model 
developed  by  Mulliken  and  Boyce  (2006).  The  Mulliken-Boyce  model  assumes  that  the 
resistance  to  the  deformation  in  thermoplastics  is  caused  due  to  intermolecular  resistance  to 
chain  segment  rotation  and  an  entropic  resistance  to  chain  alignment.  The  constitutive  model 
decomposes  intermolecular  deformation  resistance  into  contributions  from  primary  ‘a’  processes 
and  secondary  ‘/T  processes.  The  a  process  is  associated  with  rotations  of  the  polymer  main- 
chain  segments  and  the  /?  process  is  associated  with  restricted  rotation  of  side  groups.  The  model 
assumes  that  the  material  response  may  be  approximated  as  simple  superpositioning  of  the  two. 
The  model  may  be  envisioned  as  a  five  component  structure  composed  of  a  pair  of  linear  elastic 
springs  and  viscoplastic  dashpots  ‘Aa’  and  ' A ff  corresponding  to  the  a  and  processes  and  a 
nonlinear  Langevin  spring  ‘F’  acting  in  parallel.  The  basic  kinematic  assumptions  are  based  on 
multiplicative  Lee-Kroner  decomposition  of  deformation  gradient  tensor  ‘F’  applied  to  all, 

FAa  =  FlFL’  FArF^FAPr  det(FX)  =  det(F-)  =  1  (2.1) 


followed  by  an  additive  split  of  the  velocity  gradient  tensor  as, 


Li  =u-  +fplp  f;-1 

Aa  Aa  Aa  Aa  Aa 


=  L\  +LP  , 

Aa  Aa 7 


L  =  Le  +  FpLpFeJ  =  Le  +  Lp  . 

A/3  A/3  A/3  A/3  A/3  A/3  A/3 


(2.2) 


The  plastic  velocity  gradient  is  then  additively  decomposed  into  a  symmetric  stretch  (DpAa,DpAp) 
and  an  unsymmetrical  rotational  (fFpa ,  Wpfj )  component.  The  later  is  equated  to  zero  by  the 
hypothesis  of  irrotational  plasticity.  The  stretch  is  further  described  using  a  coaxial  flow  rule 


Ka 


Np  = 

Ai 


r. 

Ai 


T' 

Ai 


i  =  a,/3 


(2.3) 


where  primes  denote  deviatoric  components  of  the  Cauchy  stress  T  and  yp  the  shear  strain  rate 
which  is  given  by, 


?■  =  Y  „ ,  exp 


such  that,  T;  =  I  T\  T 


> 

i 

1 

Kr  \ 

s  +a  p 

\  i  p,'^ , 

\ |/2 

T 

i  Ai 

,  s0i=  0.077 

(2.4) 


,  and  5,  =  hj 


l 


1- 


a-O/ 


yp ;  where  yp.  is  the  pre 


(1-v,) 

exponential  factor,  AG.  is  the  activation  energy,  k B  is  the  Boltzmann  constant,  T  is  the 
temperature,  ap  .  is  the  pressure  coefficient,  p  is  the  pressure,  hj  is  the  softening  slope,  and  sj  is 
the  athermal  shear  strength  representing  the  internal  structure  of  the  material  and  is  related  to 
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shear  modulus  uj  and  Poisson’s  ratio  v,  of  the  /'th  component  through  an  appropriate  hardening 

law.  The  dependence  of  internal  structure  on  shear  strain  rate  is  a  significant  conclusion 
especially  during  loading  and  unloading  of  waves  and  this  effect  cannot  be  captured  by 
commonly  used  simple  shear  strain  dependent  stress-strain  laws  (Meyers  1994).  The  stress  in  the 
non-linear  hardening  component,  the  network  “back  stress”  due  to  the  entropic  resistance  to 
molecular  alignment,  is  taken  to  be  deviatoric  and  is  computed  as, 


r 

T  —  ^ R 

D 


3 

'''chain 


XP  \ 

chain  I 


B’ 


(2.5) 


where  Xpchain 


trace(BB) 

V  3 


is  the  stretch  on  a  chain  in  the  eight  chain  network,  lU  is  the 


Langevin  function  defined  byL(/J)  =  coth/J - ,  B\,  is  the  deviatoric  part  of  the  isochoric  left 

P 

Cauchy-Green  tensor,  Bb  =  (detFB)~2B FBFp ,  \[n  is  the  limiting  chain  extensibility  and 
CR  =  NkO  is  the  rubbery  modulus  and  N  is  the  number  of  chains  per  unit  volume,  k  is  the 
Boltzmann  constant,  and  0  is  the  absolute  temperature.  The  total  stress  in  the  polymer  is  given 
as  the  tensorial  sum  of  the  a  and  [>  intermolecular  stresses  and  the  network  (back)  stress,  that  is, 
T  =  TAa  +  TAp  +  Tb  .  Here,  we  simplify  the  composite  failure  by  proposing  a  terminal  failure  model 

for  the  bulk  matrix  and  account  for  the  damage  accumulation  in  the  failure  of  the  anisotropic 
fiber  yam.  The  bulk  matrix  polymer  resin  is  assumed  to  fail  terminally  when  failure  strain  is 
reached.  This  can  be  expressed  in  terms  of  the  molecular  chain  stretch  Xphajn  as, 


Khain*  Khain,cr’  T*Tcr’  ° ,  =  0  •  (2-6) 

Here  is  critical  molecular  chain  length  extension  Tcr  is  the  melting  temperature  and  aij  is 

internal  Cauchy  stress.  Bergstrom  et  al.  (2005)  have  demonstrated  that  this  molecular  chain 
length  failure  criterion  is  far  more  consistent  than  strain  envelopes  for  polymers.  In  addition,  the 
tenninal  damage  law  automatically  takes  into  account  the  ductile  to  brittle  transition  at  high 
strain  rates. 

The  resin  is  assumed  to  be  Epon  epoxy  and  its  properties  are  taken  from  a  recent  paper  by 
Jordan  et  al.  (2008)  and  reproduced  here.  Since  adiabatic  conditions  are  usually  assumed  during 
a  typical  high  strain  rate  events,  local  rise  of  temperature  can  be  computed  without  taking 
recourse  to  solving  a  heat  transfer  problem.  We  note  that  conversion  of  dissipative  work  to  heat 
is  a  rather  complex  process  for  polymers  where  only  a  fraction  of  post  yield  work  is  converted  to 
heat  with  the  rest  being  stored  (Garg  et  al.  2008).  In  the  current  work,  we  assume  that  the 
dissipated  work  is  contributed  solely  by  the  a  and  [>  component  whereas  the  Langevin  spring  acts 
as  pure  energy  storage  component.  If  we  denote  the  specific  heat  of  epoxy  as  Cp  and  mass 

density  as  p  we  arrive  at  the  following  relation  for  rate  of  dissipated  work  (power), 
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(2.7) 


wdiss  =  T  ■  Dp  +T  Dp 

VV  ±Aa  ■  ** Aa  T  ±  Ap  '  ^Aft 

where  overhead  dot  indicates  time  derivative  and  colon  indicates  scalar  product. 

Assuming  only  a  fraction  (i  of  this  power  is  used  to  raise  the  temperature  6 ,  we  arrive  at  the 
following  expression, 

Wdiss 

0-0  — — .  (2.8) 

PC„ 

The  exact  value  of  the  scalar  [1  varies  from  0.4  to  0.6  depending  on  strain  and  strain  rates  (Garg 
et  al.  2008).  In  the  current  problem  we  assume  the  value  of  /?=0.5. 

Table  1:  Parameters  of  the  Mulliken-Boyce  model  for  EPON  Epoxy1. 


va{d,s),  vp(6,e) 

0.38 

fo,a(6^ 

2.29xlO'V1 

Y0/0,£) 

2.0  x10  V1 

< 

3.83x  10“19  J 

< 

3.32x1  O'20/ 

a  ,  a  . 

p,a 7  p,p 

0.316 

s  ,  s  a 

ss,a 7  ss,p 

0.58  MPa 

h  ,  ha 

a 7  p 

300 MPa 

cR 

14.2  MPa 

N 

2.3  m~V2 

P 

1 1 40  Kg  cm ~3 

cP 

2000 JKg-'K' 

'Jordan  et  al.  2008 


The  strain  rate  dependent  values  of  elastic  constants  are  obtained  from  the  DMA  curves.  The 
temperature  rise  predicted  from  the  model  can  be  seen  in  the  Figure  2(a)  and  (b).  It  is  clear  from 
the  figure  that  temperature  rise  can  be  significant  and  burning  of  epoxy  is  very  likely  since  the 
ignition  temperature  of  many  epoxy  derivatives  are  about  800K.  This  prediction  is  in  agreement 
with  the  experimental  findings  of  Walter  et  al.  (2009)  who  have  reported  a  distinct  burning 
epoxy  smell  after  the  experiments.  We  will  assume  terminal  damage  of  material  occurs  at  the 
ignition  point.  Interestingly,  dependence  of  temperature  rise  on  shear  strain  indicates  that  beyond 
a  certain  point,  a  runaway  softening  is  possible  in  the  material.  This  kind  of  thermal  instability 
has  been  reported  by  in  3D-OWC  by  Pankow  et  al.  (2011)  in  their  SHPB  experiments.  However, 
we  would  like  to  note  that  in  an  actual  penetration  event  macroscopic  shear  bands  which  can 
form  underneath  a  projectile  are  quickly  interrupted  by  the  penetrating  projectile  as  well  as  other 
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failure  modes  brought  about  by  rapid  through  thickness  wave  propagation  making  post  impact 
data  reduction  difficult.  The  shear  bands  may  indeed  form  at  the  edges  of  the  projectile  due  to 
strain  localization  but  macroscopic  cracks  are  likely  to  be  interrupted  by  the  much  stronger  fibers 
present  in  all  three  directions  within  a  unit  cell.  It  also  likely  that  localized  micro  bands  can 
become  the  sites  of  void  nucleation  thus  activating  other  modes  of  mechanical  failure.  Therefore, 
possibility  of  shear  bands  alone  may  not  necessarily  cause  defeat  of  an  actual  armor. 


(b) 


Figure  2:  (a)  Temperature  vs  true  shear  strain  and  (b)  shear  stress  (Cauchy)  vs  true  shear 
strain  both  at  pure  shear  strain  rate  of 2xl06  s'1. 


3.  Constitutive  modeling  of  yarns 

Next,  we  describe  the  response  of  the  fibers.  The  fibers  themselves  are  held  together  by  the  resin 
forming  a  composite  yam.  There  yarns  are  composed  of  individual  fibers  running  axially, 
embedded  in  the  surrounding  matrix  material.  For  most  ballistic  applications  involving  3D 
composites,  S2-glass  fibers  are  commonly  employed  which  can  be  treated  as  elastic  brittle  with 
negligible  plastic  zone.  Moreover,  the  distribution  of  the  fibers  inside  the  yam  is  assumed  to  be 
uniform  for  every  cross  section  which  allows  us  to  consider  them  as  transverse  and  isotropic.  We 
assume  no  twisting  of  the  yarn  which  is  a  characteristic  of  the  most  popular  3D-OWC  produced 
by  3Tex®  Inc.  Hence,  these  yarns  can  be  treated  as  unidirectional  composite  by  themselves 
embedded  in  the  larger  3D  composite.  It  is  well  known  that  unidirectional  laminates  exhibit 
strong  nonlinearities  in  shearing  response  (Hahn  and  Tsai  1973). However,  since  the  impact 
process  causes  quick  proliferation  of  defects  and  microcracks,  all  the  nonlinearity  is  assumed  to 
come  from  damage  which  will  be  modeled  in  the  next  section.  In  addition,  the  basic  assumption 
of  yam  being  a  homogenized  continuum  is  assumed  to  be  valid  at  the  yarn  level.  This 
assumption  is  similar  to  the  one  used  by  Matzenmiller  et  al.  (1995)  and  implies  that  the  defects  in 
the  composite  material  are  treated  in  the  mathematical  model  as  having  the  equivalent  effect  on 
the  elastic  properties  as  disk-like  cracks  would  exert,  if  they  are  only  oriented  either  tangential  or 
normal  to  the  fiber  direction.  The  rate  effects  on  elastic  property  is  neglected  since  the  rate 
effects  on  the  pre-yield  modulus  is  small  as  seen  from  the  slopes  of  the  stress  strain  curve  of 
epoxy  (Jordan  et  al.  2008). 

For  moderately  high  concentration,  various  closed  form  solutions  have  been  provided  for 
fiber  composites  using  the  generalized  self  consistent  for  fiber  composites  (Christensen  1990). 
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We  use  these  closed  form  expressions  for  our  current  problem  since  they  are  simpler  and  easier 
to  use  for  various  material  parameters  compared  to  variational  methods  like  Mori-Tanaka 
estimates  (Nemat-Nasser  1993).  We  denote  the  fiber  and  matrix  Young’s  modulus,  shear 
modulus  and  Poisson  ratio  is  denoted  as  Ef,vf,iif  and  Em ,  vm ,  um  ,  respectively.  These  properties 

can  be  used  to  predict  the  five  effective  elastic  properties  of  the  transversely  isotropic  yarns  - 
En  the  axial  modulus,  v12  the  axial  Poisson’s  ratio,  ju12  the  axial  shear  modulus,  K23  the  plane 

strain  bulk  modulus  and  u2i  the  transverse  shear  modulus,  where  e1  is  the  fiber  direction.  If  we 


denote  the  fiber  volume  fraction  as  c,  we  arrive  at  the  following  expressions  (Christensen  1990), 


En  =  cE  +(l-c)Em  + 


4c(\-c)(vf-vmfnm 


(1  -c) 
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(3.3) 
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and  fi23  is  given  by  the  following  quadratic  equation, 


(  \  (  \ 

A  —  +25  ^L\  +  C  =  0 

\  Em  j  \  Em  / 


(3.5) 


where  A,  B  and  C  are  themselves  functions  of  elastic  properties  of  matrix  and  fibers  and  are 
given  as, 
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C  =  3c(l-c)2  ^-1  p-  +  ^ 

A /A 


rln—  +  —  -1  K+l 


(3.8) 


—  +  Vf  +  — ^,h-j7/  c3 

.  /  l  y,n  ) 


y,n  \Vn 

where  7jf  =  3  -  4vf  and  rjm  =  3  -  4vm  . 

Assuming  Ef=240GPa,  v/=0.3,  Em=2GPa,  vm  =  0.45 ,  c  =  0.6  and  using  the  usual 
relationship  between  elastic  constants  we  arrive  at  the  values  shown  in  Table  2. 


Table  2:  Effective  property  of  anisotropic  resin  impregnated  yarns 


En 

V12 

En 

k2. 

El3 

114.83  GPa 

0.31 

2.68  GPa 

17.16  GPa 

3.77  GPa 

3.1.  Yarn  damage  mechanics 

Predicting  damage  and  progressive  failure  in  composite  materials  under  impact  is  an  active  area 
of  research.  The  morphology  of  composite  material  induces  damage  accumulation  before 
ultimate  structural  collapse.  Hence  brittle  failure  based  criterion  will  not  yield  satisfactory  results 
as  nonlinearities  induced  by  accumulation  of  damage  would  be  stepped  over.  On  the  other  hand, 
accounting  for  every  crack  and  void  nucleated  during  loading  together  with  wave  scattering  and 
fiber-matrix  interfacial  failures  through  a  numerical  code  requires  the  kind  of  resolution  which  is 
beyond  the  reach  of  current  computing  and  imaging  technology.  Therefore  purely  computational 
approaches,  e.g.,  Oskay  and  Fish  (2007)  and  Belytschko  et  al.  (2008)  which  need  tracking  of 
cracks  explicitly  with  restrictions  on  periodicity  are  impractical  for  complex,  dynamic  loading  of 
fiber  reinforced  composites. 

The  damage  accumulation  is  assumed  to  be  completely  addressed  by  the  failure  of  the  yarn. 
Accumulating  damage  can  be  addressed  through  the  framework  of  continuum  damage  mechanics 
(CDM).  The  CDM  approach  relies  on  introducing  phenomenological  ‘internal  variables’  which 
can  track  the  degradation  of  the  material  within  the  limits  of  homogenization.  These  variables 
although  don’t  have  a  direct  bearing  to  the  micromechanics  of  crack  and  void  growth  must 
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adhere  strictly  to  restrictions  imposed  by  the  laws  of  thermodynamics.  A  good  introduction  with 
many  representative  problems  can  be  found  in  the  book  by  Lemaitre  (2005). 


Due  to  integrally  woven  geometry  of  the  3D-OWC  and  localized  ballistic  impact  loading,  we 
neglect  delamination  mode  of  damage  for  the  current  model  (Greenhalgh  2009).  The  effect  of 
friction  is  also  neglected.  The  yarn  failure  is  thus  predominantly  due  to  damage  accumulation 
through  the  intralaminar  damage  modes.  Neglecting  delamination,  we  propose  three  damage 
scalars  DF ,  Dr  and  Ds  associated  with  fiber  breakage,  transverse  cracking  and  shear  damage 

respectively.  Under  these  assumptions,  we  can  write  the  Gibbs  free  energy  for  the  yarn  as 
(Lemaitre  2005), 


G  = 
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(3.9) 


Where  otj  is  the  Cauchy  stress,  e,  is  the  fiber  direction  with  transverse  symmetry  in  the  e2  -  e, 
plane,  Ex  is  the  tensile  Young’s  modulous,  E2  =  E3  is  the  transverse  Young’s  modulous  and 

G12  is  the  shear  modulous.  We  can  use  the  above  expression  of  Gibbs  free  energy  to  obtain  the 
compliance  matrix  H  for  the  lamina, 

G  =  -—  =  [H]{o}  (3.10) 

da j 

where  e,  and  al  are  Cauchy  stress  components  in  Voigt  notation.  This  leads  to, 
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(3.11) 


The  compliance  clearly  indicates  that  the  damage  variables  are  effectively  degrading  the 
stiffness  of  the  matrix  through  the  scalar  damage  parameters.  However,  it  must  be  noted  that 
damage  caused  by  tension  can  be  different  from  that  caused  by  compression.  To  keep  track  of 
this  variation,  we  rewrite  the  damage  variables  as, 


(3.12) 
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where  (x}  is  the  McCauley  operator  defined  as  (x)  =  (x  +  |x|)/2 .  In  this  way,  the  eventual 

closure  of  transverse  cracks  under  load  reversal  is  taken  into  account.  Depending  on  the  sign  of 
the  corresponding  nonnal  stresses,  a  damage  mode  can  be  either  active  or  passive.  The  model 
also  assumes  that  the  shear  damage  variable  Ds  is  not  affected  by  the  closure  effect.  It  should  be 

noted  that  shear  damage  is  caused  mainly  by  transverse  cracks  which  do  not  close  under  shear 
stresses.  The  stiffness  degradation  through  damage  evolution  starts  as  soon  as  damage  initiation 
conditions  are  met.  The  damage  variables  then  evolve  till  a  critical  failure  value  is  reached  at 
which  point  the  material  is  said  to  have  failed  and  unable  to  resist  any  loading.  In  the  current 
work,  we  assume  that  damage  evolution  rate  depends  on  the  current  state  of  stress,  strain  and 
damage  of  that  particular  mode  with  no  interactive  modes.  The  damage  initiation  function  is 
assumed  to  depend  purely  on  stress.  Mathematically  we  can  write  this  as, 


£  =  0.  fJ°)  <  fcr' 

°  {  <Pa(cr,e,Da),  fa(o)<fcr,  (3.13) 

Da(0)  =  0,  a  =  F±,T±,S 


where  a  is  the  damage  mode,  fa,  and  (pa,  is  the  damage  initiation  and  damage  evolution 
function  associated  with  the  a  mode.  The  initiation  and  evolution  functions  represent  various 
stages  of  crack  and  void  growth  and  nucleation.  In  other  words  the  evolution  of  Da  from  0  to  a 
critical  value  Da  cr ,  represents  proliferation  of  cracks  and  voids  which  continually  degrade  the 
stiffness,  Cu  =  [//]“'  from  an  initial  value  of  C°u  to  {0}  when  all  damage  modes  have  reached 

their  critical  value  and  the  material  is  assumed  to  fail  completely.  Unfortunately  obtaining  the 
initiation  and  evolution  values  from  actual  measurements  of  crack  and  void  density,  shape,  size 
etc.  is  impossible  and  hence  phenomenological  relations  are  often  postulated.  There  is  no 
unanimity  on  the  value  of  critical  damage  variable  either.  It  represents  the  point  beyond  which 
steady  macro  cracks  cause  failure.  An  indirect  way  to  obtain  damage  evolution  and  completion  is 
using  an  energy  release  approach  wherein  energy  released  after  critical  damage  has  been  released 
is  equated  with  fracture  energy  experiments.  In  other  words,  the  area  under  the  stress- 
displacement  curve  after  the  initiation  state  is  reached  is  obtained  from  experimental  fracture 
energy  experiments.  Once  a  shape  is  postulated  (which  corresponds  to  the  damage  evolution 
law),  one  can  precisely  obtain  the  final  displacement  to  failure.  Pathological  mesh  dependency  is 
eliminated  by  normalization  through  a  characteristic  length  scale  which  depends  on  the  mesh 
size.  The  critical  damage  variable  when  the  material  is  failed  is  obtained  when  damage  energy 
dissipation  reaches  the  energy  release  rate  obtained  from  fracture  mechanics  experiments.  The 
physical  basis  of  this  theory  closely  follows  the  crack  band  theory  developed  by  Bazant  and  Oh 
(1983)  for  concrete  and  is  explained  in  detail  in  a  recent  paper  by  Lapczyk  and  Hurtado  (2007). 
The  Hashin-Rotem  criterion  (Hashin  and  Rotem  1973)  is  often  used  to  initiate  the  damage  modes 
in  fiber  reinforced  composite.  According  to  the  criteria  damage  accumulation  starts  when  scalar 
functions  representing  each  failure  mode  -  fiber  tension  ( F'f ),  fiber  compression  ( F'f  ),  matrix 
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tension  (F‘t)  and  matrix  compression  ( F'm )  reach  a  critical  value  (in  the  current  case=l).  The 
expressions  used  are  given  in  Table  3. 

Table  3:  Description  of  the  Hashin-Rotem  damage  initiation  model 


Fiber  Tension  (  on  >  0  ) 


Fiber  Compression  ( ex, ,  <  0  ) 


Matrix  Tension  ( o22  >  0  ) 


Matrix  Compression(  o22  <  0  ) 


XT  :  longitudinal  tensile  strength,  xc :  longitudinal  compressive  strength,  Y 7 :  transverse 
tensile  strength,  Yc  :  transverse  compressive  strength,  SL  :  longitudinal  shear  strength,  ST  : 
longitudinal  transverse  strength,  o'  =  He ,  H  is  the  instantaneous  degraded  compliance 
matrix  and  £  is  strain  in  Voigt  notation. 


During  ballistic  loading  however,  dynamic  conditions  are  present.  However,  since  this  failure 
model  is  used  only  for  the  yarn  and  S2  glass  fibers  exhibiting  very  little  or  no  rate  sensitivity,  the 
matrix  dominated  transverse  mode  needs  modification.  It  is  well  known  that  dynamic  crack 
fracture  is  a  complicated  phenomenon  where  fracture  energy  depends  on  duration  of  loading, 
shear  and  longitudinal  wave  speeds  as  well  as  the  possibility  of  crack  branching.  For  the  present 
case,  we  assume  that  loading  takes  place  rather  rapidly  and  as  a  first  approximation  the  dynamic 
energy  release  rate  is  approximately  25%  higher  than  the  static  case  (Meyers  1994)  for  the 
matrix  dominated  transverse  modes.  We  assume  that  transverse  failure  initiation  values  are  very 
close  to  the  yield  strength  of  polymer  matrix  resin.  Noting  the  yield  strength  increase  of  resin 
material  with  increasing  strain  rate  (Mulliken  and  Boyce  2006),  we  propose  the  following  rate 
dependent  failure  initiation  model  for  the  matrix  dominated  failure  initiation  values  YT ,  Yc  and 

SL=SL0\l  +  K 

where  Y0T ,  Y0C  and  are  respectively  the  reference  transverse  tension,  transverse  compression 
and  longitudinal  damage  initiation  strengths  evaluated  at  strain  rates  e0  and  F  is  the  deviatoric 
part  of  strain  rate  tensor.  The  reference  strain  rate,  e0  is  taken  to  be  0.001  (Jordan  et  al.  2008) 
and  A  y„ ,  a  =  T,C  and  As,  are  pre-multip lying  factor  which  indicates  the  degree  of  strain  rate 
sensitivity  for  each  mode.  The  functional  form  is  similar  to  the  one  proposed  by  Boyce  et  al. 


In-  ,  F“  =  K“  1  +  A  In-  ,  a  =  T,C 


(3.14) 
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(1988)  and  is  simplified  for  the  triaxial  loading  conditions  by  assuming  the  effective  strain  rate 
which  is  defined  as, 


e  = 


(3.15) 


Since  no  experimental  data  is  available  to  precisely  obtain  the  value  of  these  pre-multipliers 
for  3D  composites,  we  use  the  value  of  1  noting  from  Jordan  et  al.  (2008)’s  uniaxial  experiments 
that  yield  strength  roughly  doubles  when  strain  rate  increases  by  two  orders  of  magnitude. 
However,  experimental  testing  can  be  done  to  yield  a  more  accurate  value  of  this  exponent.  The 
failure  initiation  values  used  in  this  work  are  given  in  Table  4. 


Table  4:  Failure  initiation  values  for  anisotropic  yarns. 


XT  (MPa) 

Xc  (MPa) 

Y0t  (MPa) 

Y0C  (MPa) 

S'(]  (MPa) 

6999 

8999 

6999 

399999 

7999 

Finding  the  compression  fracture  energy  for  the  model  remains  a  challenge  because  not  many 
standardized  tests  have  been  done  for  various  fiber-matrix  combinations.  The  only  tests  available 
in  literature  which  can  be  directly  used  for  our  damage  model  have  been  done  by  Pinho  et  al. 
(2006)  using  a  compact  tension  and  compact  compression  specimen  for  T300/Epoxy  composite. 
We  have  adapted  Pinho  and  coworker's  fracture  data  and  the  experiments  done  by  Parhizgar  et 
al.  (1982),  using  values  from  0°  and  90°  in  the  fracture  energy  calculations  of  Pinho  et  al.  (2006). 
Using  these  approximations  we  arrive  at  the  following  damage  properties  for  the  yarn  (Table  5): 


Table  5:  Fracture  energy  released  rate  for  various  damage  modes 


Gft,c  (N/mm) 

Gfc  ,c(N/mm) 

Gmt,c(N/mm) 

Gmc,c  (N/mm) 

(Fiber-tension) 

(Fiber-compression) 

(Matrix-tension) 

(Matrix-compression) 

91699 

76999 

26899 

22299 

4.  Constitutive  modeling  of  transversely  isotropic  yarns 

Fibers  themselves  are  held  together  by  the  resin  forming  a  composite  yam.  There  yarns  are 
composed  of  individual  fibers  running  axially,  embedded  in  the  surrounding  matrix  material.  For 
most  ballistic  applications  involving  3D  composites,  S2-glass  fibers  are  commonly  employed 
which  can  be  treated  as  elastic  brittle  with  negligible  plastic  zone.  Moreover,  the  distribution  of 
the  fibers  inside  the  yarn  is  assumed  to  be  uniform  for  every  cross  section  which  allows  us  to 
consider  them  as  transverse  and  isotropic.  We  assume  no  twisting  of  the  yarn  which  is  a 
characteristic  of  the  most  popular  3D-OWC  produced  by  3Tex®  Inc.  Hence,  these  yarns  can  be 
treated  as  unidirectional  composite  by  themselves  embedded  in  the  larger  3D  composite. 

The  body  starts  from  initial  configuration  Bo  at  temperature  9o  and  at  the  end  of  loading 
(lasting  for  time  At)  is  at  the  new  configuration  B  with  the  new  temperature  9.  We  instead 
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envision  the  step  through  a  couple  of  intermediate  stress  free  states  denoted  by  B*  and  B.  There 
is  no  heat  transfer  to  the  surroundings  from  the  body. 


Figure  3:  Decomposition  of  deformation  through  intermediate  configuration. 


Figure  4:  Intermediate  decomposition  on  temperature-specific  entropy  diagram.  Reversible  processes  are 
shown  in  solid  arrows  and  irreversible  processes  in  dashed  lines. 

The  intermediate  states  can  be  described  in  a  thennodynamic  specific  entropy  (entropy  per 
unit  reference  volume)-temperature  (r|-9)  diagram  as  show  in  Figure  2  below.  From  the  picture  it 
is  clear  that  entropy  is  produced  during  the  plastic  step  (B0-B*)  due  to  various  inelastic  processes 
like  motion  of  polymer  chains,  crack  fonnation  etc.  We  consider  this  entropy  to  be  much  smaller 
in  comparison  to  other  subsequent  entropic  processes.  We  consider  the  next  thermal  step  (B*-B) 
to  be  reversible  (equilibrium)  and  purely  thennal  taking  the  temperature  from  initial  Go  to  9i.  The 
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heat  transferred  is  supplied  by  the  plastic  dissipation  from  the  previous  plastic  step  (B0-B*).  The 
next  step  (B-B)  can  be  combination  of  reversible  and  irreversible  step  depending  upon  the 
presence  of  viscous  effects. 


4.1.  Thermodynamics  of  intermediate  steps 

In  this  section,  we  look  in  detail,  the  thermodynamics  of  the  intennediate  steps. 

Step  Bo-B*:  This  step  is  an  isothermal,  isochoric  and  adiabatic  elastic  unloading.  The  work  in  this 
step  is  used  to  create  inelastic  phenomena  like  dislocations,  polymer  chain  sliding,  cracking  etc. 
The  total  work  done  is  composed  of  elastic  part  and  an  inelastic  part.  The  elastic  part  is  assumed 
to  increment  (or  decrease)  the  internal  energy  whereas  part  of  the  plastic  part  is  assumed  to 
produce  heat. 

Thus  thermodynamically  this  can  be  written  as: 


Ae  ,  =  -dw  ,  +8q  .  +  3q 

B„-B  B„-B  1B„-B  1B.-B 


(4.1) 


where  Ae  denotes  internal  energy  change  per  unit  initial  (reference)  volume,  dw  denotes 
specific  work  done  by  the  system,  dq  denotes  specific  thermal  heat  added  to/generated  by  the 
system,  other  irreversible  sources  such  as  acoustic/light  etc  are  denoted  by  dq  at  this  point. 
The  total  work  during  this  elastic  unloading  can  be  divided  into  elastic  and  plastic  work.  The 
elastic  work  goes  to  change  the  internal  energy  and  part  of  the  plastic  work  causes  heat 
production  and  other  irreversible  loss.  Thus: 


Ae 


Bn-B 


=  ~dwt-B'  +Sq 


Bn-B 


lBn-B 


lB„-B 


(4.2) 


If  we  assume  that  a  certain  fraction  k  of  plastic  work  goes  directly  to  heat  we  have: 

SC/b0-B'  =K(dWB0-B')  =  ,C{<5wPl) 


(4.3) 


where  5wpl  is  the  plastic  work. 


Step  B*-B:  In  this  step,  the  heat  produced  in  the  previous  step  causes  internal  energy  change  and 
no  stress  is  applied  to  the  system.  Thus,  it  is  a  workless  step.  Reapplying  the  first  law  of 
thennodynamics  once  again,  we  have: 

AtV-«  =~dWB'-B  +dC<B-B  +  ^P-B  (4'4) 


Now  dw  ,  =  dq  ,  =0  due  to  purely  thermal,  stress  free-state  assumption.  Assuming  that  the 

total  thennal  heat  in  step  B0-B*  is  applied  during  step  B*-B,  we  have  dO  ,  =  dO 

^B  —B  ^~^Bq  —B 

Dividing  throughout  by  reference  volume  V  =  V . ,  we  get  dq  ,  =  dq  , .  Thus,  we  have: 

Bq  B  1 B  —B  1  Bq  —B 

=  SCl,f  -B  =  SClB0-B'  =  K(dwP‘  )  (4‘5) 
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From  basic  thermodynamics,  we  have  at  constant  stress: 

q  .  -  C„0  =  kwp!  (4.6) 

1B0-B  S  v  ' 

where  Cs  is  the  specific  heat  per  unit  reference  volume  and  wpI  is  the  plastic  power.  Thus 
temperature  equation  can  be  written  as: 

f'Cs  dd  =  KwplAt  (4.7) 

J  60 

where  A t  is  the  time  increment.  If  specific  heat  per  unit  current  volume  is  used,  the  thennal 
volume  change  must  be  taken  into  account. 

4.2.  Computing  stress  from  elastic  step 

The  elastic  step  consists  of  the  process  denoted  by  B-B.  This  can  be  reversible  in  the  case  of 
perfectly  elastic  material  (which  has  no  entropy  sources)  or  partly  irreversible  in  case  of 
viscoelastic  materials.  This  step  can  also  be  purely  isothermal  if  elastic  and  viscous  heating  are 
completely  neglected.  We  will  proceed  with  each  step  relaxing  the  conditions  on  the  way  to 
derive  the  most  general  formulation.  However,  the  thermomechanical  building  blocks  are  built 
here  which  would  be  applicable  later. 

We  follow  the  standard  Coleman-NoU  formulation  of  extracting  thermo-mechanical 
governing  equations  from  the  intrinsic  dissipation  inequality: 

Dini=^S:Ce  -e  +  9r]>  0  (4.8) 

where  S  is  second  Piola-Kirchoff  stress  and  Cf  is  the  elastic  part  of  left  Cauchy-Green  stress,  e 
and  rj  are  specific  internal  energy  and  specific  entropy  per  unit  referential  volume,  respectively, 
and  6  is  the  temperature.  If  the  process  is  purely  reversible  the  equality  holds  or  else  the 
inequality. 

It  is  clear  from  the  above  fonnulation  that  derivatives  of  energy  and  entropy  must  be 
calculated  with  respect  to  time.  These  time  derivative  are  often  converted  to  strain  derivatives 
through  chain  rule.  We  carry  out  systematic  analysis  of  these  derivatives  below. 

We  multiplicatively  decompose  the  elastic  defonnation  gradient  as  following: 

Fe=jfFe,  det(Fe)  =  l  (4.9) 

Similarly,  the  left  Cauchy-Green  stress  can  be  written  as: 

Ce  =  FeJF\  Ce  =  J2e/3Ce,  det (Ce)  =  1  (4.10) 

This  multiplicative  split  in  the  defonnation  gradient  can  be  augmented  by  an  additive  split  in 
the  specific  internal  energy  and  specific  entropy.  For  equilibrium  processes  (no  viscoelasticity), 
this  can  be  written  as: 
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e(0.C')-eM(0,J.)+et.(0.C)  (4.11) 

and 

n(0.c')-ti«(e.j.)+va,(e,c)  (4.12) 

The  energy  derivatives  can  be  written  as  following: 


p*  = 


to*  \ 


S.  =2 


(  de.  \ 

1  iso  1 

dCe 


(4.13) 


It  must  be  noted  that  p,  is  not  pressure  (that  would  be  derivative  with  respect  to  Helmholtz  free 
energy  at  constant  temperature).  Similarly  S*  is  not  the  usual  isochoric  stress. 


The  computation  of  entropic  derivative  is  eased  by  the  thermodynamic  quantity  called  the 
Gruneisen  tensor.  The  Gruneisen  tensor  is  a  thermodynamic  measure  of  variation  of  entropy  with 
strain  and  is  a  useful  quantity  in  further  derivation.  This  tensor’s  decomposition  into  volumetric 
and  isochoric  (deviatoric)  parts  will  greatly  aid  in  computation  of  time  derivative  of  entropy 
which  is  an  essential  quantity  for  future  derivations.  The  Gruneisen  tensor  is  written  as: 


r  = 


d)j  \ 
dC 


drivol{(hJe)\  .(dr,iso(d,Ce)\ 


dCe 


dCe 


It  can  be  shown  that  this  tensor  can  be  written  as: 


(4.14) 


r  = 


c. 


y  —  ce~l  +  j_2/3r :  P7 
2 


(4.15) 


d%, 


^Viso _ J  r»r  t  ^ 


where  y  =  —7^,  T  =  and  P7  =  I  ~^Ce  1  ®Ce  is  the  transpose  of  the  projection  tensor. 


dJ. 


8Ce 


4.2.1.  Purely  elastic  without  structural  heating 

We  are  now  in  a  position  to  compute  the  stress  expressions.  We  start  with  the  simplest  case  of 
isothermal  purely-elastic  step  with  no  heating.  In  this  case,  the  temperature  in  the  step  (B-B) 
remains  constant  at  0\.  The  intrinsic  dissipation  for  this  step  can  be  written  as: 

D,„-^:C--e+e,n  (4.16) 


This  can  be  written  using  chain  rule  as: 


-S:&- 

de  .  ■  de. 

VOl  J  |  ISO 

\ 

■.te 

+  0, 

(Ho,  T  ,  drhso.^e\ 

2 

dJ  e  dCe 

\  e 

J 

1 

dJ  e  dCe  ' 

\  e  / 

(4.17) 
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This  can  be  recast  as: 


D  =-S:Ce- 

mt  2 


/ 


pje  +S*  1  ~ 


&\e. ' 


Ce  \ 

yj  +r :  — - 
"  2 


V 


(4.18) 


-^S:C'-(p,-Y8l)jr-Us.:Cr  +  m':C- 


From  tensor  analysis  we  know: 

re 

J  =J  CeA  :  — 

e  e  2 

_2 

te  =  2J~^T  :  — 
e  2 


(4.19) 

(4.20) 


Thus  we  get  the  expression  for  intrinsic  dissipation  as: 

D,„  -  fs :  C-  -  (p.  -  re,  )JC :  £1  -  jj  (s.  -  flf) :  Tr : 


2 


ce 

T 


)  r 


(4.21) 


Now  due  to  perfect  elasticity  with  no  heating,  the  dissipation  must  disappear,  leading  to: 

_2  _2 

S=(p.-y0 \  )jeCe~ 1  -  Jp  (s.  -  6X  f) :  Pr  =  (p  -  yd,  )jeCe~l  -  Jp  P  :  (s.  -  6,  f )  (4.22) 


The  conversion  to  Cauchy  stress  can  be  done  using  the  following  expression: 
o  =  —FSFeT 

Je 

or 


(4.23) 


o  =  — 
J 


z 

(p.  -j/0,  )jeFeCe-lFeT  -JpFeV-.( S.  -^f)Fer 


Recalling  the  deviatoric-isochoric  split  of  the  defonnation  gradient,  we  arrive  at: 


o-(p.-ya 1  )/+</;1Fe[P:(s.-01r) 


(4.24) 


(4.25) 


4.2.2.  Purely  elastic  with  structural  heating 

In  the  thermo-elastic  step,  there  is  temperature  change.  Thus  the  time  rate  of  change  of  internal 
energy  and  entropy  can  be  written  as: 
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e(6,Ce)  = 


dj 


\ 


j  + 


( 


h 


de,JdjCe) 

dCe 


\ 


/e 


:Ce+  —  6 

1  dd  Le 


and 


(4.26) 


j  + 

e 

e 

t>nje,C)\ 

\  dJe  ) 

{  da  ) 

la#) 

6  '  f 

mce)= 


Thus,  the  intrinsic  dissipation  would  become: 


(4.27) 


D  =-S:Ce- 

mt  2 


( 


In  other  words: 


1 


/ 


/  \  ^ 

l  dc  i 

—  o 

/ 

+e 

\dd  )  e 

\  /c  y 

\ 

yj  +r  :  — +  (—  I  e 

2  30  L 


(4.28) 


D  =  —S :  Ce  - 

mt  2 


P*J  e  +S*  1  Y 


re  \  ( 


+e 


ce  \ 
yj  +T:  — 
e  2 


t 


e 


I  drj  j  |  de  ' 

\dd j  \d0y 


0 


(4.29) 


Clearly,  the  last  tenn  which  is  coefficient  of  temperature  rate  is  identically  zero  from  basic 
thennodynamics.  Thus,  the  form  of  equation  developed  in  the  earlier  section  is  preserved: 


D  =-S:Ce- 

mt  2 


t 


-  Ce^  ( 


ptJe+K-  — 


+e 


re  \ 
yj  +T:  — 
2 


(4.30) 


This  is  similar  to  previous  section  and  thus  yields  the  following  Cauchy  stress  expression: 


o  =  {P,-yd)l  +  J-elFc\  P:(S.-6T) 


?eT 


(4.31) 


However,  an  important  distinction  is  that  the  temperature  is  no  longer  constant  and  its  evolution 
with  time  must  be  separately  provided  for  by  a  temperature  evolution  equation. 

The  temperature  rise  in  this  step  is  now  a  combination  of  both  structural  elastic  as  well  as 
inelastic  effects.  The  temperature  evolution  can  be  derived  as  follows: 


r](Ce,d)  =  -^-:&  +  ^-9 

dce  dd 

We  recall  the  following  thermodynamic  identities: 

2  (  dp 


(4.32) 


r  = 


c 


dC 


Cc=o[-\ 

\dd)c 


(4.33) 


Thus  the  above  equation  can  be  written  as: 
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(4.34) 


C  C 

i)(Ce,6)  =  ^T:Ce  +  ^6 


e 


This  can  be  written  as: 


1 


6i)  =  -6CcX-.Ce  +  Cce6 

The  heat  transfer  equation  can  next  be  written  as: 
61 7  =  -V  • q  +  D  +R 


(4.35) 


(4.36) 


where  q  is  the  heat  flux,  Dint  is  the  internal  dissipation  and  R  is  sum  of  local  volumetric  heat 
sources.  Since  the  process  is  perfectly  elastic  and  adiabatic,  Dint=  0,  R= 0,  q=().  Thus  eliminating 
entropy  among  the  equations,  we  get: 


0  =  ^dCceT:Ce  +  Cced 

Thus,  temperature  evolution  is  given  by: 

e  =  --6T:Ce 
2 

It  can  be  shown  that: 

0  =  -9T  :  FeTDeFe 


(4.37) 


(4.38) 


(4.39) 


where  De  +  F1 )  is  the  elastic  stretch  rate  given  in  terms  of  the  velocity  gradient  V  . 


If  the  temperature  at  the  end  of  this  process  is  d2  the  Cauchy-stress  would  be: 


0=(p*-Yd2)i+J:xFe\v-.(  S.-02r) 


?eT 


(4.40) 


4.2.3.  Visco-elastic  with  structural  heating/cooling 

Viscoelastic  phenomena  is  an  inelastic  phenomena  which  can  occur  even  without  a  threshold 
effect.  One  way  to  model  it  is  to  assume  viscoelasticity  to  come  from  purely  non-equilibrium 
(NE)  processes  which  is  additive  to  the  steady  state  internal  energy  and  entropy  through  internal 
variables  as: 


e{d,CeA  )  =  e”(6,Ce)  +  evisc(0,?l 


(4.41) 

(4.42) 
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where  e°°  is  the  steady  state  response  whereas  evlsc  is  the  non  equilibrium  response  which  is 
assumed  to  depend  only  on  the  deviatoric  part  of  strain.  Here,  we  would  need  an  additional 
evolution  equation  for  the  internal  variable  is: 


i  =  a,p  =  \...N 


Writing  the  dissipation  inequality: 


1 


Din,  =  —S :  Ce  -  e(d  ,Ce ,...,U  +  di)(d,Ce 


(4.43) 


(4.44) 


or 


where 


Dini=l-S:Ce -e»{9,C)  +  evisc{d,Ce  ,...,%N)  +  d[i)°°(d,Ce)  +  i)''isc(d,Ce,l  ,...,^)]>0  (4.45) 


„  c)pvisc  .  r)Pvisc 


se 


Let  us  assume  a  simple  evolution  law: 


(4.46) 


£  +  —  =  /3  Ce 

~  a  •  a 

X 


(4.47) 


plugging  in  above  yields: 


de 


vise  {  .  &  \ 


similarly: 


dL 


p  ce 


a  ) 


d£. 


dev,sc  . 

+ - :d 

86 


n  vise 

+?2—.e 


a  \ 


J  a  ) 


Thus,  the  dissipation  inequality  becomes: 

Dmt  =  \s-.a-  e-(d,ce)+dine,ce)-2Pa 


dd 


d(evisc-6r rn.^e 


d£ 

n 


,  d(evisc  -  6rjvisc)  rg 


dL 


:  —  >  0 
T 


(4.48) 


(4.49) 


(4.50) 


Clearly  the  viscous  terms  multiplying  the  deviatoric  strain  rates  have  units  of  stress.  In  addition, 
recognizing  the  form  of  Helmholtz  free  energy,  we  get: 
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(4.51) 


1  s)iliVlSC  & 

D.nt  =  -  S :  Ce  -  e*  (6,  Ce )  +  &rf  ( 6,Ce)~  Svisc  :—  +  ^  ~7—  ■  —  s  0 

2  2  0  x" 

a  a 


where 


^  /  VWC  _  vise  \ 

j~-2 


d^’ 


(4.52) 


3#a  ^  ^ 

assuming  from  construction  of  appropriate  function  the  following  is  made  to  identically  be 

di/jvisc  p 

satisfied  ^  ^ — :  —  >  0  and  noting  from  equation  (4.20),  we  get: 


31 o  *„ 


D  =  — 

m'  2 


5-J_35v“c:Pr 


) 


:Ce-r(0,Ce)  +  0?)"(0,Ce)^O 


(4.53) 


Now  the  above  can  be  further  simplified  as: 


1  „  .  4 

D.nt  =  ^S:Ce  -e”(0,Ce)  +  0?f(0,Ce);>O,  S  =  S-Je3Svbc:  VT 


(4.54) 


This  enormously  simplifies  our  analysis  for  we  can  carry  out  our  usual  derivations  for  the 
remaining  part  thereby  arriving  at: 


_2  _2 

S  =  (P:-  yd)Je  C'-'  -  J  3  (s:  -  6t ) :  PT  =  ( pt  -  yd)Je  Ce~l  -  J3P :  (SC  -  Of) 


(4.55) 


where  p™  = 


ax 


Soo  _ 

* 


dC 


.  The  conversion  to  Cauchy  stress  can  be  done  using 


equation  (4.23)  as  follows: 


a-  —  Fc 

J 


/  -2 

J  3SV,SC :  Pr 


FeT  ={p.  -yd)l  +  J;'Fe  [P:(s.  -6r) 


(4.56) 


or 


a  =  (pt.-yd)l  +  J-lFe\P:(st-dT  +  Svisc ) 


?eT 


(4.57) 


The  temperature  rise  in  this  step  is  now  a  combination  of  both  structural  elastic  as  well  as 
inelastic  effects.  The  temperature  evolution  can  be  derived  as  follows: 


i)  (d,Ce,Zl,...,ZN)  =  -^-:Ce  +  ^-:6  +  yN 

1  N  dCe  86  Z'«-'  drE  a 

a 


(4.58) 
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We  recall  the  following 
equation  can  be  written  as: 


r  = 


c 


d/J 

lie7 


}e 


drj\  dij 

dd)  C  d£a 


thus  the  above 


drf 


dL 


■L 


(4.59) 


This  can  be  written  as: 


dr)  =  -dC  T:Ce  +C  ed  +  dyN  :£ 

'  2  c  c  £>a-\  " 

The  heat  transfer  equation  can  next  be  written  as: 
61 7  =  -V  •  q  +  D.nt  +  f? 


(4.60) 


(4.61) 


where  q  is  the  heat  flux,  D,nt  is  the  internal  dissipation  and  R  is  sum  of  local  volumetric  heat 
sources.  Since  the  process  is  adiabatic,  q=0.  Thus  eliminating  entropy  among  the  equations,  we 
get: 


D  +R  =  idC  T:Ce  +C  ed  +  &yN  :£ 

mt  2  C  C  1  Qp  a 

-a 

Thus,  temperature  evolution  is  given  by: 

C  e9  =  D  +R-idC  eT:Ce -ey" 

Ce  I"'  2  C  Z>a= 1  Qh  « 

Plugging  in  the  value  of  intrinsic  dissipation: 


C£  =  eJ 

C  Aj  a= 1 


n  dr/' 


n 


„VISC  \  i 

:£  +R--BC  .r  :C’-eJ" 

£ja= 1  Qg  a  9  Ce  Zja= 1  '  r- 


\n  drf 


a  / 


n 


or 


c  ed  =  R--6C  J:Ce~yN  S  :| 

Ce  O  Ce  Z^a=l  « 


(4.62) 


(4.63) 


(4.64) 


The  temperature  rise  can  be  computed  if  the  evolution  equations  for  the  internal  variable  are 
known: 


The  final  Cauchy  stress  then  can  be  computed  at  final  temperature  02  as: 


a  =  (p,-y62)l  +  J-elF‘ 


P:|St-02r 


ieT 


(4.65) 


(4.66) 
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4.3.  Fiber  composite  -  hyperelastic  case 

Let  us  explore  the  case  of  fiber  composite  in  detail.  Let  us  denote  the  initial  orientation  of  the 
fiber  by  the  unit  vector  a0 .  Thus  the  internal  energy  can  be  written  as: 


e(0, Ce )  =  e{0, u0,Ce)  =  e(6, /,  ,/2,/3,/4,/5)  =  e(d,  Je ,I{  ,I2 ,14 ,15) 
where  Ik  and  Ik  are  the  strain  invariants  of  the  given  problem  and  are  given  by: 

/,  =  tr(Ce),I2  =^(lf-tr(Ce2)),I2  =det  (Ce),I4  =a0-Cea0,I5  =  a0-Ce2a0 
and 

l  -tr(C-),l  -i(7,2-(r(C-2)),74  -ae-C-a0J,  =«o'C'2«„ 

Assuming  a  deviatoric/isochoric  split  in  the  energy: 

e(d,Je,h  J2  JA  Js  )  =  evol(0’Je)+eiso(0Jl  Jl  Ja  JS  ) 

where 

evol(6,Je)  =  K(6)(Je- 1)2 

The  fictitious  modified  isochoric  stress  S,  can  be  written  as: 


g  O  deiso  2\N  9eiso  .  91  a 

/ j  a=l,< 


dC  l.a-3  dC 

=  YxI  +  Y2C  +  Y4do  ®d0  +  y5  (d0  ®Cd0  +  Cd0  ®du ) 


where  y{  =  2 


1 5e  -  de ......  \ 


Kdl, 


ISO  _|_  J  ISC 


dl 


de:. 


de,.„ 


de:. 


V2  =-2-=^,  Ya  =2 -^,y5  =2—= 


2  ) 


dl 


dl. 


dl 


In  indicial  notation, 


X  d0  ®  Cd0 ,  Xtj  a0iCjka0k  aoiaok^jk 
X  —  d0C  ®d0,  Xy  =  a0kCkia0 }  =  a0ka0jCki 


Thus: 


(4.67) 


(4.68) 


(4.69) 


(4.70) 


(4.71) 


(4.72) 


(4.73) 


^*ij  ~  Yi  dy  +  y2Cij  +  y4 a0ja0 j  +y5  {a0ia0kCjk  +  a0ka()jCki  j  (4.74) 

The  anisotropic  potential  takes  the  following  form  (simple  extension  of  Mooney-Rivlin): 

eiSO{dJ{  ,I2J4,l)  =  aimi  -3 )  +  ff2(60(72  -3)  +  a4(d){J4  -l)  +  a5(0)(75  -1)  (4.75) 

Therefore, 
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(4.76) 


ft  =2 [ax{d)+Ixa2(d)),  y2  =-2 a2{d\  y4  =2 a4(0),  y5  =2 a5(0) 
For  a  simple  isotropic  case: 
l  de, 


p«  = 


~'vol 


\  dJe  )e 


=  2K(6){Je  - 1) 


S.  =  2 (al(0)  +  Ila2(0))l- 2 a2 (0)C  =  2a 1  ( 0)1  +  2a2 (6>)(/,  I - C) 
For  pure  translation  F  =  /.  C  =  C  =  I ,  /,  =3 


S.  =  ft  I  +  y2I  +  y4a0  ®a0  +  y5  (d0®a0+a0®a0) 

=  (ft  +ft2  )I+(n  +2  ft  )«0®«o 
=  KXI  +  K2a0  ®a0 
=  KJ  +  K2A0 


Hence, 


P  :  S.  =  Kx  xO  +  ftx 


A- 


=  X 


4o- 


(4.77) 

(4.78) 


(4.79) 


(4.80) 


4.4.  Fiber  composite  -  viscoelastic  case 

The  next  case  study  pertains  to  composites  with  finite  viscoelasticity.  The  internal  energy 
function  is  assumed  to  be  of  the  following  fonn: 


e(e,c%l . ^N)  =  e(6,d0,Ce,^ . £n) 


(4.81) 


For  the  sake  of  brevity,  let  us  assume  that  is  only  one  internal  variable  for  the  time  being,  thus: 


e(0,Ce,£l,..,£N)  =  e(0,do,Ce,i D 


(4.82) 


From  the  discussion  of  the  previous  section,  it  is  transparent  that  all  we  need  is  a  way  to  estimate 
the  viscous  stress  through  the  internal  variable.  The  dissipative  Helmholtz  potential  can  be 
written  as  xjjv,sc{6 ,£)  and  the  evolution  equation  is  given  by: 


Thereby  the  viscous  stress  is  given  by: 


dljj' 


d£a 


=  2 fi 


dip' 


(4.83) 


(4.84) 
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A  simple  form  of  free  energy  that  satisfies  the  second  law  dissipation  inequality  Sv,sc :  %  >  0  is 
given  by: 


V* (4.85) 

where  Cvisc  is  a  viscosity  tensor.  An  easier  simplification  would  be  if  it  is  assumed  to  be  isotropic 
thereby  giving: 

(4-86) 

Thus: 

Svisc=C°v,J  (4.87) 


For  this  simple  case,  one  can  write  the  evolution  equation  directly  in  terms  of  viscous  stress  as: 


q  vise 


(4.88) 


4.5.  Yarn  damage  mechanics 

Predicting  damage  and  progressive  failure  in  composite  materials  under  impact  is  an  active  area 
of  research.  The  morphology  of  composite  material  induces  damage  accumulation  before 
ultimate  structural  collapse.  Hence  brittle  failure  based  criterion  will  not  yield  satisfactory  results 
as  nonlinearities  induced  by  accumulation  of  damage  would  be  stepped  over.  On  the  other  hand, 
accounting  for  every  crack  and  void  nucleated  during  loading  together  with  wave  scattering  and 
fiber-matrix  interfacial  failures  through  a  numerical  code  requires  the  kind  of  resolution  which  is 
beyond  the  reach  of  current  computing  and  imaging  technology.  Therefore  purely  computational 
approaches,  e.g.,  Oskay  and  Fish  (2007)  and  Belytschko  et  al.  (2008)  which  need  tracking  of 
cracks  explicitly  with  restrictions  on  periodicity  are  impractical  for  complex,  dynamic  loading  of 
fiber  reinforced  composites. 

Continuing  from  the  previous  grant  period,  we  now  extend  the  continuum  damage  mechanics 
framework  to  include  damage  mechanisms  suitable  for  the  high  strain  rate  regime. 

The  total  internal  energy  without  damage  can  be  written  as: 


e(0,  Je ,  1 1  ,I2,l 


4  5 


h  )  =  F 


vol 


{d,Je)  +  eiso(d,Ix,I2,h,I5) 


(4.89) 


We  assume  that  the  volumetric  and  isochoric  damages  are  a  result  of  different  mechanisms,  with 
volumetric  part  directly  contributed  by  pressure  and  the  isochoric  part  through  shear.  Thus,  the 
energy  functional  can  be  written  as: 


e  =  (l -D  ,)e  A6,J  )  +  e  (6,1,  Dl  ,D 2  ,D3  ,D4  ) 

V  vol  A  vol  \  5  eJ  iso\  ’  1  ’  2  ’  4  ’  5  ’  iso 5  iso 5  iso  ’  iso  A 


(4.90) 
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where  Dvoi  and  Diso  are  the  volumetric  and  shear  damage  parameters  respectively.  We  postulate 
the  volumetric  damage  evolution  law  as  follows: 


D  =DN)lN 2 

vol  vol  cr 


(  „  \ 


V  °  I 

\  ecl  / 


(4.91) 


where  N\,  Ni  and  W,  are  empirical  exponents,  lcr  is  the  empirical  characteristic  flaw  size,  p  is 
pressure  and  oeq  is  the  von-Mises  stress  at  the  material  point.  The  isochoric  damage  variables 
D1  ,  D2  D  ,  D4  represent  uniaxial  tension  (fiber  parallel),  compression  (fiber 

perpendicular),  transverse  shear,  and  in-plane  shear,  respectively  and  each  have  their  own 
evolution  rules. 


5.  Numerical  Examples 

There  have  been  relatively  few  experiments  conducted  to  quantify  ballistic  penetration  data  for 
3D  orthogonal  woven  fiber  composites.  One  of  the  preliminary  experiments  on  the  system 
currently  being  studied  has  been  done  by  Gama  et  al.  (2001).  The  authors  carried  out  bullet 
impact  experiments  on  3D-OWC.  The  panels  were  obtained  from  3Tex®  Inc.  In  their 
experiments,  bullets  of  5  mm  diameter  made  from  Steel  were  fired  orthogonally  at  these  panels. 
The  authors  reported  photographic  evidence  as  well  as  measurements  of  incident  and  exit  speeds. 
In  this  section  we  will  simulate  the  event  and  discuss  the  results  obtained  from  the  model  against 
experimental  results. 

5.1.  Composite  plate  geometry 

The  particular  three-dimensional  woven  fiber  composite  geometry  considered  here  in  constructed 
by  repeating  a  fundamental  unit  cell,  which  is  shown  in  Figure  5(a)-(c)  and  is  developed  after 
Chou.  The  overall  dimensions  of  the  unit  cell  are  9.23 6mm  x  9.932 mm  x  5.390 mm.  The  fill  (X), 
warp  (Y)  and  Z  tows  are  idealized  as  prismatic  bodies.  The  fiber  tows  are  discretized  using  finite 
strain  continuum  shell  elements  and  the  resin  by  fully  integrated  three  dimensional  hexahedral 
elements.  A  minimum  of  four  continuum  shell  elements  are  used  to  discretize  the  smallest 
dimension  of  a  tow.  In  the  overall  discretization,  the  aspect  ratios  of  the  elements  are  kept 
similar.  A  total  of  127616  of  these  continuum  shell  elements  were  used  each  with  thickness 
0.1mm  and  5  integration  points  used  for  through  thickness  integration  to  model  the  X,  Y  and  Z 
yams.  A  total  of  47424  of  3D  fully  integrated  hexahederal  elements  were  used  to  discretize  the 
resin  and  projectile.  This  mesh  density  has  been  kept  the  same  for  all  analysis  performed  in  this 
paper. 

The  plate  is  clamped  at  the  boundaries  allowing  for  no  slip  at  the  clamps.  Furthermore, 
symmetry  boundary  conditions  are  used  to  model  only  quarter  of  a  plate  and  hence  4  repeating 
unit  cells  as  our  computational  domain  can  be  used  to  model  a  plate  4  times  as  large. 
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Figure  5:  (a)  Impact  system  setup  with  a  quarter  plate  and  quarter  of  a  rigid  right  cylindrical  projectile  (b) 
unit  cell  (c)  meso  scale  components  making  up  the  unit  cell. 

Since  we  are  interested  in  investigating  penetration,  a  sufficiently  high  enough  projectile 
speed  will  be  taken.  Penetration  is  a  complex  process  and  depends  on  a  large  number  of 
parameters.  Projectile  momentum  as  well  as  mechanical  properties  detennine  the  time  of 
penetration.  At  high  speed  penetration  boundary  conditions  are  of  little  significance  due  to  highly 
localized  damage.  In  this  section  we  present  some  numerical  examples  of  simulations  of 
projectile  impact  at  500  m/s  on  a  three  dimensional  woven  fiber  composite  plate.  The  projectile 
is  assumed  to  be  a  rigid  right  circular  cylinder  under  isothermal  conditions  with  mass  of  1  gram 
and  diameter  of  5 mm  to  represent  the  projectile  used  in  the  ballistic  experiment. 

5.2.  Experimental  validation  and  key  insights 

5.2.1.  Hyperbolic  damage  profile 

The  damaged  zone  which  corresponds  to  the  region  of  high  stress  transients  has  a  hyperbolic 
front  centered  at  the  point  of  impact.  Clearly  the  evolving  stress  contours  in  Figure  6(c)  indicate 
an  advancing  hyperbolic  stress  front.  Our  simulations  from  Figure  6  clearly  reproduce  the 
characteristic  diamond  shaped  failure  region  at  the  front  and  at  the  back  as  observed  in  the 
experiments.  Our  simulations  also  show  localized  Z-fiber  damage  as  reported  by  the  experiments 
(Gama  2001). 


Figure  6:  (a)  Experimental  photographs  of  plate  rear  side  by  Gama  et  al.  (2001);  (b)  qualitative  comparison 
of  simulation  results  with  experiments  done  by  Gama  et  al.  (2001)  (c)  hyperbolic  profile  of  the  advancing 
stress  front. 

5.2.2.  Highly  Directional  Rear  Face  Tensile  Damage  Distribution 

The  next  figure  shows  the  spatial  distribution  of  the  tensile  damage  field  of  the  rear  face  just 
after  a  perforation  event  has  taken  place.  In  the  experimental  results,  there  is  a  marked  difference 
in  the  damage  directionality  of  the  strike  face  compared  to  the  rear  face  as  seen  in  Figure  7(a). 
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The  Figure  7(b)  corroborates  the  fact.  The  directionality  is  defined  with  respect  to  the  direction 
of  Z  fiber  windings  over  the  matrix. 


(b) 


Figure  7:  (a)  Tensile  damage  region  just  after  penetration  (b)  fiber  damage  and  pullout  visible  in  experiments 

done  by  Gama  et  al.  (2001). 

5.2.3.  Bulk  Matrix  Rate  Effects  Show  Two  Distinct  Regimes 

Disappearance  of  rate  effects  (Figure  8(a)-(b))  is  very  interesting  as,  theoretically,  rate  effects 
should  be  minimal  for  quasi-static  indentation  and  maximum  for  high  speed  impacts.  We  explain 
this  apparent  paradox  by  realizing  that  as  the  matrix  resin  stiffens  at  high  strain  rates,  there  is  a 
competing  nonlinearity  due  to  the  onset  of  damage.  Hence,  corresponding  to  a  high  speed 
impact,  the  improvement  in  ballistic  resistance  due  to  rate  effects  is  effectively  nullified  by  a 
rapid  onset  of  material  failure.  This  anomaly  of  diminishing  rate  effects  have  been  observed  by 
before  during  the  shock  plate  experiments  of  (Pankow  et  al.  2011)  and  cannot  be  explained  by 
hot-spot  analysis  only  without  incorporating  a  damage  law  into  the  formulation. 


Figure  8  (a)  Projectile  velocity-time  profile  comparison  between  models  with  and  without  rate  effects  at 
500/nA  (b)  Projectile  velocity  comparison  between  models  with  and  without  rate  effects  at  100/n/s 

5.2.4.  Damage  Anisotropy  Effect  of  Z-fiber  on  Matrix  and  Fiber  Damage 

The  effect  of  Z-fibers  is  more  than  simply  preventing  delamination  across  inter-ply  layers.  Figure 
9(a)  shows  that  the  Z-fibers  running  on  the  surface  of  the  Y-fibers  in  an  orthogonal  direction 
scatter  the  stress  waves  away  from  them,  thereby  reducing  the  immediate  area  of  damage.  The 
damage  is  thus  restricted  primarily  in  the  Y-direction.  Z-fibers  thus  cause  directional  damage  in 
the  fibers  but  diffuse  damage  in  the  matrix  in  the  yarns.  Figure  9(a),  (c)  and  (b),  (d)  show  the 
fiber  compression  and  tensile  damage  damage  in  3D-OWC  and  2D  layered  composites, 
respectively.  The  left  and  right  panels  of  Figure  10(a)-(f)  show  that  in  spite  of  highly  localized 
tenninal  the  matrix  damage  at  Z  fiber  sites  has  potential  for  void  growth  and  interlaminar  cracks 
thereby  reducing  the  ability  of  Z-fibers  to  act  effectively  during  future  impact  events.  Clearly  the 
purported  benefits  of  Z  fibers  over  stitched  composites  are  that  it  doesn’t  physically  damage  the 
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composite  structure.  In  deed  this  conclusion  has  been  experimentally  observed  by  Pankow  et  al. 
(2011)  during  shock  loading  of  these  panels  who  noted  that  beyond  6%  Z-fiber  concentration, 
the  stiffness  properties  appreciably  deteriorate. 


Figure  9:  (a):  Strike  face  fiber  tension  damage  for  3D  OWC  (b)  Strike  face  fiber  tension  damage  for  layered 
composite  (c)  Strike  face  fiber  compression  damage  for  3D  OWC  (d)  Strike  face  fiber  compression  damage 
for  layered  composite 


Figure  10:  (a):  Strike  face  fiber  tension  damage  for  3D  OWC  (b)  Strike  face  fiber  tension  damage  for  layered 
composite  (c)  Strike  face  fiber  compression  damage  for  3D  OWC  (d)  Strike  face  fiber  compression  damage 
for  layered  composite 

Finally,  Table  6  gives  the  summary  of  the  challenges  and  achievements  for  the  macro-scale 
modeling. 


Table  6:  Summary  of  challenges  and  their  resolution 


Modeling  and  Computational  Challenges 

Damage  laws 

Damage  initiation  rate  sensitivity 
Resin  ductile  to  brittle  transition 
Damage  evolution 


Successful  Resolution _ 

Meso  level  micromechanics 
Modified  Hashin 

Terminal  thermo-mechanical  damage  law 
Crack  band  type  regularized  fracture  mechanics 
based  model 


5.3.  Study  of  anisotropic  composite  material  response  with  changing  material  parameter 

We  now  examine  the  effect  of  Gruneisen  tensor  on  the  composite  materials  response.  The 
Gruneisen  parameter  is  an  important  material  property  quantifying  the  entropy  dependence  on 
deformation  for  a  material.  In  a  typical  shock  physics  application  (Davison  2008),  this  is  treated 
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as  a  scalar  variable.  This  is  justified  through  the  reasoning  that  shock  regime  comprises  of  bulk 
response  of  a  material  which  is  well  past  the  shear  limit  of  a  typical  material.  In  our 
computational  material  model,  we  found  that  for  anisotropic  materials,  the  full  tensorial 
representation  is  critical  to  capture  the  essential  physics  of  the  process.  In  the  Figure  1 1  and  12, 
we  show  that  varying  the  anisotropy  of  Gruneisen  tensor  can  have  significant  impact  on  the 
incipient  plasticity  of  a  material.  We  show  that  as  this  tensor  is  rotated  about  [100]  (Figure  11) 
and  [001]  (Figure  12),  the  J?  changes  significantly  even  though  the  bulk  response  may  not.  We 
also  report  that  even  in  the  case  of  isotropic  solids  with  anisotropic  Gruneisen  tensor,  same 
effects  show  up. 


(a)  (b)  xlO4 


(c) 


o3 

PL, 


S 


a 


(d) 


Figure  11:  Anisotropic  material  response  with  changing  Gruneisen  tensor  (rotation  about  [  100]). 
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(C)  (d) 


Figure  12:  Anisotropic  material  response  with  changing  Gruneisen  tensor  (rotation  about  [001]). 


6.  Summary 

We  achieved  or  even  exceeded  key  goals  to  achieve  a  multiscale  damage  model  for  the  impact 
problem  proposed.  The  meso  scale  damage  mechanics  avoids  the  theoretical  pitfalls  associated 
with  regular  multiscale  methods  based  on  fiber  level  classical  homogenization  as  well  as  unit  cell 
level  averaging  methods.  The  anisotropic  yam  properties  are  modeled  using  generalized  self- 
consistent  method  for  fiber  composites.  This  method  has  the  advantage  of  yielding  closed  form 
expressions.  The  model  predicted  various  experimental  observations  from  direct  penetration  tests 
as  well  as  standardized  split  Hopkinson  and  shock  plate  impact  tests  as  well.  We  found  that  both 
Z-fiber  and  bulk  matrix  rate  effects  can  improve  ballistic  resistance  depending  on  the  momentum 
of  the  projectile  and  dimension  of  the  plate.  However,  in  the  penetration  range,  competing 
nonlinearities  tend  to  nullify  each  other  producing  both  localized  effects  and  potential  crack 
nucleation  sites  through  wave  propagation.  We  also  found  that  the  effect  of  Z-fiber  on  damage 
distribution  is  not  straightforward  and  analysis  of  both  liber  and  matrix  damage  mode  explained 
many  of  the  unique  damage  features  of  3D-OWC.  Although  obtaining  data  for  fitting  damage 
models  remains,  we  have  drastically  cut  the  number  of  arbitrary  parameters  necessary  to  obtain 
predictive  capabilities.  This  model  of  the  yam  is  further  extended  into  a  large  deformation 
model.  We  propose  a  damage  law  based  on  consistent  thermomechanics  for  a  typical 
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unidirectional  composite,  which  can  be  easily  extended  for  composites  with  greater  degree  of 
anisotropy.  We  examined  the  effect  of  material  parameter  such  as  Gruneisen  tensor  on  the 
composite  materials  response.  In  our  computational  material  model,  we  found  that  for 
anisotropic  materials,  the  full  tensorial  representation  is  critical  to  capture  the  essential  physics  of 
the  process.  We  also  report  that  even  in  the  case  of  isotropic  solids  with  anisotropic  Gruneisen 
tensor,  same  effects  show  up.  However,  we  could  not  validate  our  model  due  to  lack  to 
experimental  data.  Overall,  this  mesoscale  model  while  less  computationally  expensive  than  a 
standard  homogenization,  it  still  remains  computationally  prohibitive  for  large  plate  impacts.  The 
next  step  in  addition  to  parameter  estimation  would  be  to  develop  a  continuum  level  model  to 
make  the  model  amenable  to  large-scale  computations. 
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